#########################################################################
# Time-stamp: <05/08/17 11:51:56 ostolop>
#########################################################################
package EP::EPC::Core::Statistics;

use strict;

use EP::Config;
use EP::Common;
use EP::Common::General;
use EP::Common::Error ':try';

use EP::EPC::Core::Clustering::exp;

use EP::EPC::Core;
use EP::EPC::Core::Component;

use EP::Rclient;

use File::Temp;
use File::Basename;

our @ISA = qw ( EP::EPC::Core::Component );

sub process_component {
  my $self = shift;

  my ( $pval_plot_url, $pval_table_url ) = $self->run_t_test()
    if $self->{query}->param('process_component_id') == 104;
  my $process_result = {
    pval_plot => {
      type  => "image/png",
      value => $pval_plot_url,
      label => "15 lowest p-values plot"
    },
    pval_table => {
      type  => "xml/table",
      value => $pval_table_url,
      label => "Genes and their p-values"
    },
    src_dataset_id => $self->{src_dataset}->id,
    dst_dataset_id => $self->{src_dataset}->id
  };

  $self->{CD}->{process_result} = $process_result;
  return $self->record_analysis_history();
}

sub run_t_test {
  my $self = shift;

  throw EP::Common::Error( -value => AFFY_TRANSFORM_RAW_DATA )
    if $self->{src_dataset}->type eq "affy";

  my $ttest_param;

  my $dataset_folder_path =
    $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $self->{src_dataset}->folder->location;
  my $bin_file   = $dataset_folder_path . "/" . $self->{src_dataset}->filename . ".bin";
  my $outputfile = new File::Temp(
    DIR    => $dataset_folder_path,
    SUFFIX => ".ttestpvals",
    UNLINK => 0
  );

  # one class
  if ( $self->{query}->param('class_setup') eq 'one_class' ) {
    my $control_genes    = $self->{query}->param('control_genes');
    my $control_gene_ids =
      EP::EPC::Core::Clustering::exp::get_gene_ids( $self->{src_dataset}, $control_genes,
      $self->{query}->param("rows_regexp") )
      if $control_genes;
    my $control_mean = $self->{query}->param('control_mean');

    $$control_gene_ids[$_]++ foreach 0 .. $#{$control_gene_ids};

    $ttest_param = {
      html_root      => $EP::Config::h->{HTML_ROOT},
      src            => $bin_file,
      filename       => $self->{src_dataset}->filename,
      outputdir      => $dataset_folder_path,
      outputfile     => "$outputfile",
      controlgeneids => defined $control_gene_ids
      ? join( ",", @$control_gene_ids )
      : '',
      controlmean => $control_mean,
      alpha       => $self->{query}->param('alpha'),
      correction  => $self->{query}->param('correction'),
    };
  }

  # two classes, one file
  elsif ( $self->{query}->param('class_setup') eq 'two_class_one_dataset' ) {
    my @class_one = parse_col_ranges( $self->{query}->param('class_one') );
    my @class_two = parse_col_ranges( $self->{query}->param('class_two') );
    $ttest_param = {
      html_root      => $EP::Config::h->{HTML_ROOT},
      src            => $bin_file,
      filename       => $self->{src_dataset}->filename,
      outputdir      => $dataset_folder_path,
      outputfile     => "$outputfile",
      classone       => join( ",", @class_one ),
      classtwo       => join( ",", @class_two ),
      controlgeneids => '',
      controlmean    => 0,
      alpha          => $self->{query}->param('alpha'),
      correction     => $self->{query}->param('correction'),
    };
  }
  elsif ( $self->{query}->param('class_setup') eq 'two_class_two_datasets' ) {
    my @pnames      = $self->{query}->param;
    my @dataset_ids =
      map { /chk_dataset_(\d+)/ } grep { /^chk_dataset/ } @pnames;
    my @action_ids =
      map { /chk_action_(\d+)/ } grep { /^chk_action/ } @pnames;

    die "Can't specify more than two datasets for t-test!"
      if scalar @dataset_ids + scalar @action_ids > 2;

    my ( $dataset1, $dataset2 );
    if ( scalar @dataset_ids == 2 ) {
      $dataset1 = EP::Common::DBI::Dataset->retrieve( $dataset_ids[0] );
      $dataset2 = EP::Common::DBI::Dataset->retrieve( $dataset_ids[1] );
    }
    elsif ( scalar @dataset_ids == 1 and scalar @action_ids == 0 ) {
      $dataset1 = $self->{src_dataset};
      $dataset2 = EP::Common::DBI::Dataset->retrieve( $dataset_ids[1] );
    }
    elsif ( scalar @dataset_ids == 1 and scalar @action_ids == 1 ) {
      $dataset1 = EP::Common::DBI::Dataset->retrieve( $dataset_ids[0] );
      $dataset2 = EP::Common::DBI::Action->retrieve( $action_ids[0] )->dst_dataset;
    }
    elsif ( scalar @dataset_ids == 0 and scalar @action_ids == 1 ) {
      $dataset1 = $self->{src_dataset};
      $dataset2 = EP::Common::DBI::Action->retrieve( $action_ids[0] )->dst_dataset;
    }
    elsif ( scalar @action_ids == 2 ) {
      $dataset1 = EP::Common::DBI::Action->retrieve( $action_ids[0] )->dst_dataset;
      $dataset2 = EP::Common::DBI::Action->retrieve( $action_ids[1] )->dst_dataset;
    }

    my $dataset1_binfile =
        $EP::Config::h->{EPNG_USRDATA_PATH} . "/"
      . $dataset1->folder->location . "/"
      . $dataset1->filename . ".bin";
    my $dataset2_binfile =
        $EP::Config::h->{EPNG_USRDATA_PATH} . "/"
      . $dataset2->folder->location . "/"
      . $dataset2->filename . ".bin";

    $ttest_param = {
      html_root      => $EP::Config::h->{HTML_ROOT},
      src            => $dataset1_binfile,
      src2           => $dataset2_binfile,
      filename       => $self->{src_dataset}->filename,
      outputdir      => $dataset_folder_path,
      outputfile     => "$outputfile",
      classone       => '',
      classtwo       => '',
      controlgeneids => '',
      controlmean    => 0,
      alpha          => $self->{query}->param('alpha'),
      correction     => $self->{query}->param('correction'),
    };
  }

  EP::Rclient->R(<<__RCODE__);
app_root = "$EP::Config::h->{ROOT}";
source("$EP::Config::h->{ROOT}/cgi-lib/tmpl/R/ep.run_t_test.tmpl.R");
ep.ttest( "$ttest_param->{src}"
        , "$ttest_param->{src2}"
        , c($ttest_param->{classone})
        , c($ttest_param->{classtwo})
        , $ttest_param->{alpha}
        , "$ttest_param->{correction}"
        , c($ttest_param->{controlgeneids})
        , $ttest_param->{controlmean}
        , "$ttest_param->{outputfile}" );
__RCODE__

  my $pvalsfile       = $outputfile;
  my $output_www_base =
    "$EP::Config::h->{EPNG_USRDATA_URL}/" . $self->{src_dataset}->folder->location . "/";

  my $pvals_errbar_file_url = $output_www_base . basename($outputfile) . "_errbar.png";

  my $pvals_table_file = new File::Temp(
    DIR    => $dataset_folder_path,
    SUFFIX => ".ttestpvals.xml",
    UNLINK => 0
  );
  my $pvals_table_file_url = $output_www_base . basename($pvals_table_file);

  my $xml = qq{<?xml version='1.0'?>};

  if ( -s $pvalsfile ) {
    open P, $pvalsfile or die "Can't open $pvalsfile for reading";

    $xml .= "<table>";

    $xml .= "<head>";
    $xml .= "<data>$_</data>"
      foreach ( "Gene", "(adjusted) p-value", "mean/difference of means", "confidence interval" );
    $xml .= "</head>";

    while (<P>) {
      chomp;
      my @d = split /\t/;
      @d = ( $d[0], $d[1], $d[3], "($d[2], $d[4])" );

      $xml .= "<row>";
      $xml .= "<data>$_</data>" foreach @d;
      $xml .= "</row>";
    }
    close P;
  }
  else {
    $xml .=
        "<table><head><data>No genes found with p &lt; "
      . $self->{query}->param('alpha')
      . ".</data></head>";
  }

  $xml .= "</table>";

  print $pvals_table_file $xml;

  return ( $pvals_errbar_file_url, $pvals_table_file_url );
}

# Sub: parse_col_ranges
#   Parses column ranges, meant to take comma-separated strings like 1,2,3-40,51 and unparse these
sub parse_col_ranges {
  my $in = shift;

  $in =~ s/-/\.\./g;    # change 1,2,3-40,51 into 1,2,3..40,51 (perl array specification)

  my %x;
  @x{ eval($in) } =
    undef;   # now eval $in to make it into array and make it into keys of a hash, to remove repeats

  return sort { $a <=> $b } keys %x;
}

1;

